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FIELD OF THE INVENTION 

The present invention relates to a method for forming more rapidly a stochastic 
numerical model of Gaussian or related type, representative of a porous heterogeneous 
medium (such as a hydrocarbon reservoir for example) calibrated in relation to data 
5 referred to as dynamic, characteristic of the displacement of fluids in the medium such 
as, for example, production data (pressures obtained from well tests, flow rates, etc.). 

The method according to the invention finds applications in the sphere of 
underground zone modelling intended to generate representations showing how a 
certain physical quantity is distributed in a zone of the subsoil (permeability notably), 
10 best compatible with observed or measured data, in order for example to favour the 
development thereof. 

BACKGROUND OF THE INVENTION 

Optimization in a stochastic context consists in determining realizations of a 
stochastic model which meet a set of data observed in the field. In reservoir 

15 engineering, the realizations to be identified correspond to representations, in the 
reservoir field, of the distribution of carrying properties such as the permeability or 
porosity. These realizations form numerical reservoir models. The available data are, for 
example, isolated permeability or porosity measurements, a spatial variability model 
determined according to isolated measurements or data directly related to the fluid 

20 flows in an underground reservoir, i.e. pressures, breakthrough times, flow rates, etc. 
The latter are often not linearly related to the physical properties to be modelled. A 
randomly drawn realization is generally not in accordance with the whole of the data 
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collected. Coherence in relation to the data is integrated in the model by means of an 
inverse procedure : 

- Tarantola, A., "Inverse problem theory - Methods for data fitting and model parameter 
estimation", Elsevier Science Publishers, 1987. 

5 The simplest technique is therefore the trial-and-error method. This approach 

consists in randomly taking realizations until a realization meeting the data is obtained. 
It affords the advantage of conservation of the spatial variability model but it requires a 
prohibitive calculating time. It is therefore very rarely used in practice. 

An option that is often preferred is based on the gradients calculation. The gradients 
10 methods allow to modify an initial realization in a direction of search which is estimated 
from the gradients. The modifications are applied iteratively until data calibration is 
considered to be acceptable. The gradients methods are attractive because of their 
efficiency. However, they are no longer suitable when the number of parameters, i.e. the 
number of permeability and porosity values forming the numerical model, is large. 
15 Besides, they do not allow to modify the realizations while respecting the spatial 
structure of the stochastic model. 

More recently, a geostatistical parameterization technique has been introduced to 
constrain, by gradual deformation, the stochastic realizations to data on which they 
depend non-linearly. It is described in patents FR-2,780,798 and FR-2,795,841 filed by 
20 the applicant, and in the following publications, notably : 

- Hu, L.Y., 2000, Gradual deformation and iterative calibration of Gaussian-related 
stochastic models : Math. Geology, Vol.32, No.l, 
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- Le Ravalec, M. Et al., 2000, The FFT moving average (FFT-MA) generator : An 
efficient numercial method for generating and conditioning Gaussian simulations : 
Math. Geology, Vol.32, No.6, 

- Hu, L.Y., Blanc, G. And Noetinger, B. (2001) : Gradual deformation and iterative 
5 calibration of sequential stochastic simulations. Math. Geology, Vol.33, No.4. 

This method has been successfully applied in various cases, notably from data 
obtained from oil development fields, as described in the following documents : 

- Roggero, F. et al., 1998, Gradual deformation of continuous geostatistical models for 
history matching, paper SPE 49004 : Proc. SPE Annual Technical Conference and 

10 Exhibition, New Orleans, 

- Hu, L.Y. et al., 1998, Constraining a reservoir facies model to dynamic data using a 
gradual deformation method, paper B-01 : Proc. 6 th European Conference on 
Mathematics of Oil Recovery (ECMOR VI), 8-1 1 September 1998, Peebles, Scotland. 

As reminded in detail hereafter, the gradual deformation method allows to gradually 
15 modify a realization of a stochastic model of Gaussian or Gaussian-related type while 
respecting the spatial structure of the model. 

Stochastic optimization 

Let yob. = (^ot* f*» 3 .„ t y M <*») be the field data and / = (fu /*•••. /m) the 
corresponding responses simulated for a realization z of a given stochastic model Z. In 
20 general, the responses / = (fu U) are obtained by solving numerically the direct 
problem. Thus, if z represents a permeability field, data / can be pressure 
measurements. In this case, they are simulated from a flow simulator. The goal of a 
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stochastic optimization is to produce realizations of Z which reduce the differences 
between the observed data and the numerically simulated corresponding responses. 
These differences are measured by the following objective function : 

^ m=l 

5 Coefficients co m are weights assigned to data W*". f m are functions of realization z. 

In this sense, minimization of the objective function is a problem with several variables. 

Let N be the number of grid cells or of components forming realization z. N is often 
very large (10 4 ~ 10 7 ). It is therefore very difficult to perform an optimization directly in 
relation to the components of z. Furthermore, realization z, even modified, must remain 
10 a realization of Z. The gradual deformation method allows these difficulties to be 
overcome. 

Random search from the gradual deformation method 

We now consider a random field Z that can be transformed into a Gaussian random 
field Y. The gradual deformation technique allows to construct a continuous chain of 
15 realizations by combining an initial realization y 0 of Y with another realization u,, 
referred to as complementary, of Y, u, being independent of y 0 (Figure la). The 
combination coefficients are for example cos(t) and sin(t), and the combined realization 
meets the relation : 

y(t) = y 0 cos t + Ui sin t 
20 where t is the deformation parameter. 

Another realization chain construction technique consists in combining the initial 
realization no longer with one, but with P complementary realizations Up (p=l,P) of Y 
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(Figure lb). The coefficients of the combination are such that the sum of their squares is 
1. 

Once the chain is formed, it can be explored by varying deformation parameter t 
and one tries to identify, from among all the realizations of this chain, the realization 

5 which minimizes the objective function. This minimization is performed in relation to t. 
Parameterization according to the gradual deformation method allows to reduce the 
number of dimensions of the problem from N to 1, where N is the number of values that 
constitute the field to be constrained. Furthermore, the sum of the combination 
coefficients squared being 1, the optimized realization yCW) still is a realization of Y : it 

10 follows the same spatial variability model as all the realizations of Y. 

However, if the exploration of the realizations space is restricted to a single chain, 
our possibilities of sufficiently reducing the objective function are greatly limited. The 
above procedure therefore has to be repeated, but with new realization chains. These 
realization chains are constructed successively by combining an initial realization which 
15 is here the optimum realization determined at the previous iteration, with a 
complementary realization of Y, randomly drawn each time. Thus, at iteration /, the 
continuous realization chain is written as follows : 
y, (t) = y M cos t + u, sin t . 

y„ is the optimum realization defined at iteration /-l and the u, are independent 
realizations of Y. The latter are also independent of y 0 . This formulation implies that the 
realization chain corresponds to a hyperellipse of dimension N. 

Minimizing the objective function in relation to t allows to improve or at least to 
preserve calibration of the data each time a new realization chain is explored. This 
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iterative minimum search procedure is continued as long as data calibration is not 
satisfactory. The random side of the method lies in the fact that, upon each iteration, the 
complementary realization is drawn at random. In fact, the direction of search that is 
followed from the optimized realization at the previous stage is random. The direction 
of search, for a given chain and from the optimum realization defined above, is : 



dy,(t) 
dt 



- -_)/,_, sin 0 + u, cos 0 
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This direction of search only depends on u,. Furthermore, u, being independent of 
the complementary realizations already generated Ui, u 2 ,..., u M and also of y 0 , the 
direction of search at the start of each new chain is orthogonal to the tangent defined for 
the previous chain at the same point (Figure 2). Although it may seem appropriate to 
select a direction of search orthogonal to this tangent, there is an infinity of possible 
orthogonal directions. 

Experience shows that, after several iterations, the new directions of search no 
longer contribute significantly to the decrease of the objective function (Figure 6). 

It has also been considered combining the initial realization not only with one, but 
with several complementary realizations. In this case, the number of deformation 
parameters increases : it is equal to the number of complementary realizations involved 
in a gradual combination. Although the optimization process is then more flexible, 
several parameters have to be managed, which is not easy. Besides, such a process is 
not necessarily more efficient because it can depend on the execution of a larger number 
of direct flow simulations. 
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SUMMARY OF THE INVENTION 

The method according to the invention allows more rapid formation of a numerical 
model representative of the distribution of a physical quantity in a porous heterogeneous 
medium such as an underground zone (oil reservoirs, aquifers, etc.), constrained by data 
5 collected in the medium (dynamic data characteristic of the displacement of fluids in the 
medium), collected by measurements (in production, injection or observation wells for 
example) or previous observations. 

It comprises an iterative process of gradual deformation wherein an initial 
realization of at least part of the selected model of the medium is linearly combined 

10 with at least a second realization independent of the initial realization, the coefficients 
of this linear combination being such that the sum of their squares is 1, and an objective 
function measuring the difference between a set of non-linear data deduced from said 
combination by means of a medium simulator and said geologic and dynamic data is 
minimized by adjusting the coefficients of the combination, the iterative process being 

15 repeated until an optimum realization of the stochastic model is obtained. 

The method is characterized in that the rate of gradual deformation to the optimum 
model representative of the medium is accelerated by selecting as the second realization 
to be combined with the initial realization at least one composite realization obtained by 
selecting beforehand a direction of descent defined as a function of the gradients of the 
20 objective function in relation to all the components of said initial realization. 

The composite realization is obtained for example by linear combination of a set of 
independent realizations of the model, the coefficients of the combination being 
calculated so that the direction of descent from the initial realization y is as close" as 
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possible to the one defined by the gradients of the objective function in relation to all 
the components of the initial realization. 

Optimization is for example carried out from a deformation parameter which 
controls the combination between the initial realization and the composite realization. 

In cases where said combination affects only part of the initial realization, the 
iterative process of gradual deformation is applied to a Gaussian white noise used to 
generate a Gaussian realization and the derivatives of the objective function with 
respect to the components of the Gaussian white noise are determined. 

According to an implementation mode, the initial realization is combined with a 
certain number M of composite realizations, all obtained by composition from P m 
independent realizations of Y, the optimization involving M parameters. 

In other words, the method essentially comprises a new gradual combination 
scheme which takes account of the information provided by the gradients at the initial 
point of any chain of realizations. Construction of a chain always starts from an initial 
realization and from a set of complementary realizations, all independent and coming 
from the same stochastic model. The initial realization is however not directly combined 
with the complementary realizations. The complementary realizations make it possible 
to explore the realizations space in different directions. These directions are not 
equivalent : some allow to get nearer to the optimum. At this stage, a realization 
referred to as composite realization is elaborated by combining the complementary 
realizations only. A chain of realizations is then created from the initial realization and 
from this composite realization. This chain, like the chain proposed in the basic case of 
gradual deformation, can be explored by means of a single deformation parameter. 



The composite realization is constructed so as to propose a direction of search along 
the chain as close as possible to the direction of descent given by the gradients. As 
mentioned above, all complementary reactions are not equivalent : the composite 
realization takes the best of each complementary realization. 

The method allows to reach the formation of the numerical model representative of 
the medium more rapidly. 

BRIEF DESCRIPTION OF THE FIGURES 

Other features and advantages of the method according to the invention will be 
clear from reading the description hereafter of an application given by way of non 
limitative example, with reference to the accompanying drawings wherein : 

- Figures la, lb show gradual deformation schemes (referred to as GDM) that are 
already known, 

- Figure lc shows the gradual deformation scheme (GBC) corresponding to the method 
according to the invention, 

- Figure 2 shows realization chains in a Euclidean space with N dimensions, where the 
tangent at the level of the optimized realization for a realization chain 7-1 (RC M ) is 
orthogonal to the direction of search for the initial realization of the next realization 
chain / (RQ), 

- Figure 3 shows the projection of the direction of descent v in the subspace U defined 
by the orthonormal base formed by P independent realizations (u,,..., Up), 
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- Figures 4 A to 4C show permeability distribution examples respectively for the 
reference realization, the initial realization and the realization constrained to the 
pressure data, 

- Figures 5A to 5E respectively show the variations as a function of time of the 
bottomhole pressures respectively simulated for the five wells (BHP-OBS1, BHP- 
OBS2, BHP-PROl, BHP-OBS3, BHP-OBS4) of Figures 4, respectively for the 
reference (data), initial and constrained (match) permeability distributions, and 

- Figure 6 shows the evolution (OF) of the objective function as a function of the 
number k of flow simulations performed, GDM1 corresponding to an optimization 
carried out by combining an initial realization and a single complementary realization, 
GDMGBC3, to an optimization carried out by combining the initial realization and a 
composite realization constructed from three complementary realizations, GDMGBC10, 
to the optimization carried out by combining the initial realization and a composite 
realization constructed from ten complementary realizations, and GDMGBC30, to an 
optimization carried out by combining the initial realization and a composite realization 
constructed from thirty complementary realizations. 

DETAILED DESCRIPTION 

The method according to the invention allows to orient upon each iteration the 
construction of the realization chain in order to reach a desirable direction of search. 
The technique selected takes advantage of the information provided by the direction of 
descent defined by the gradients of the objective function J. It can be implemented by 
means of a numerical simulator of a type known to the man skilled in the art, such as the 
ATHOS or ECLIPSE simulators. 
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Search oriented by the direction of descent 

We consider at this stage the direction of descent (evolution towards a minimum 
value, whether local or not) defined by the gradients of the objective function J in 
relation to all the components of realization z. These gradients are deduced from 
5 sensitivity coefficients df m I dz, : 

3J / r robs \ Qfm 

The problem of sensitivity coefficients calculation has been widely dealt with in the 
scientific literature. The following document can for example be referred to : 
- Sun, N.Z., Inverse problems in groundwater modelling, Kluwer Acad. Publ., 
10 Dordrecht, The Netherlands, 1994. 

In particular, the adjoint state technique allows to calculate all of these coefficients 
by solving the direct problem and its adjoint problem, within a time interval equivalent 
to twice the time required for solution of the direct problem. 

Since the goal is to integrate the information provided by the gradients in the 
15 gradual deformation method, we have to go back to the Gaussian realization which 
results from the transformation of z. Two cases can be distinguished according to 
whether isolated conditioning data, i.e. measurements of z at certain points, are 
available or not. 

When no data are available, if y is the Gaussian realization obtained by 
20 transforming z, then : 

dJ = dJ dz n 
6y„ dz„dy n 
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In the opposite case, we assume that z is known at certain points a. Thus, the 
Gaussian realization s, deduced from the transformation of z, is a conditional realization 
of Y. To obtain s, we generate a non-conditional realization y of Y and condition it to 
the values known at a by kriging. The conditional realization is deduced from : 

s = s* + (y-y*). 

s* and y* are respectively calculated by kriging from the real data and the data of y 
generated at the level of points a. We can then show that : 

dJ dz. 



dJ 
dy„ 



ds n 

0, Vw = a 



For continuous physical properties, dVdy,, or dz^dSn are calculated from the 
10 anamorphosis function allowing to transform realization z into a Gaussian realization. 
When the physical properties considered are category or discrete properties, these 
derivatives do not exist. The gradients techniques then cannot be applied. 

These various relations are of direct interest if a realization just has to be deformed 
globally. On the other hand, if it is desired to deform it locally, the gradual deformation 
15 method has to be applied to the Gaussian white noise x used to generate y. In this case, 
an additonal stage is required : calculation of the derivatives of the objective function 
with respect to the components of the Gaussian white noise. 

To illustrate the calculation of these derivatives, we suggest to concentrate on the 
particular case where the non-conditional Gaussian realization y is obtained from the 
20 FFT-MA generator described in the article published by Le Ravalec et al. 2000 
mentioned above. 
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The basic principle of this FFT-MA (FFT-Moving Average) generator is to 
transform a Gaussian white noise x into a Gaussian realization y correlated from a 
convolution product : 



y = g*x. 



Function g results from the expansion of covariance function C such that C g * 
where g* is the transpose of g. The derivatives of the objective function with respect to 
the components of the Gaussian white noise are : 

a/_ _Y pJ dy, 

The discrete expression for the convolution product leads to y, = 2>,-*** , which 

10 implies ay/at = g*. If we introduce this formulation in the derivatives of the objective 
function, we show that : 

dx n fdy, {ty )„ 

This formula expresses the fact that the derivative of the objective function with 
respect to the n* component of the Gaussian white noise is given by the n* component 
15 of the field obtained by convolving all the derivatives of the objective function with 
respect to the components of the Gaussian realization with the kernel of the covariance 
function. From the framework established for the FFT-MA generator, these derivatives 
are determined as follows. 

1-Calculation of the Fourier transform of dJ / dy, these derivatives being obtained 
20 by means of the direct numerical simulator ; 
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2- Multiplication of this Fourier transform with that of g that is provided by FFT- 
MA; 

3- Calculation of the inverse Fourier transform of the previous product. 

The time required for calculation of these derivatives is negligible it represents an 
additional time of 2/3 in relation to the simulation of a Gaussian realization by FFT- 
MA 

Whatever the realization generator and the direct numerical simulator, we assume 
hereafter that we can define a direction of descent from the derivatives of the objective 
function. If optimization of the objective function is performed in relation to this 
direction of descent only, the coherence of the realization in relation to the spatial 
variability model is generally destroyed. In the section hereunder, we integrate the 
information provided by the derivatives of the objective function in the gradual 
combination scheme. 

Taking account of the derivatives of the objective function in the gradual 
15 deformation process 

We consider the realization chain y,(t) constructed from y 0 and from another 
realization u of Y (Figure la). Now, instead of using a complementary realization u as it 
is, we randomly draw P complementary realizations Up (p=l,2,...,P) of Y and we 
substitute u for a combination of the iw (Figure lc). This combination is not any 
20 combination : it follows the following construction mode : 
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The resulting realization u is a realization of Y and it is independent of y 0 . 
According to Eq.l, u is also the direction of search calculated for chain y,(t) at the 
starting point y 0 . We therefore suggest to construct u so that direction u is as close as 
possible to the direction of descent given by the gradients at y 0 . 

We first define the space V of the vectors with N dimensions provided with the 
scalar product : 

y„ and yi „ are respectively the n* components of vectors y s and yj . Let U be a 
subspace of V defined by the orthonormal base (u,^,...^). The direction of search in U 
10 which is the closest to the direction of descent v is given by the projection of v in U 
(Figure 3) : 

v y =Z( v ' M *W 

p 

By normalizing this vector, we obtain the desired direction u. The combination 
coefficients X of Eq.(2) are thus : 



The realization u thus defined is referred to as composite realization. 

We have considered so far the construction of realization chains by combination of 
the initial realization with a composite realization. This technique can however be 
generalized to the construction of chains involving a certain number M of composite 
20 realizations, all obtained by composition from P m realizations of Y, these ££4 P. 
realizations being independent. Optimization then involves M parameters. This 
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technique increases the degree of freedom in the optimization process, but it requires 
management of M optimization parameters. 

Numerical example 

We construct a synthetic reservoir model on which the method according to the 
invention is tested. 

The synthetic reference reservoir is shown in Figure 4A. It is a monolayer reservoir 
comprising 51x51 grid cells which arelO m thick and 40 m in side. The permeability 
distribution is lognormal with an average of 200 mD and a standard deviation of 100 
mD. The logarithm of the permeability field is characterized by an isotropic spherical 
variogram and a correlation length of 480 m. The other petrophysical properties are 
constant : the porosity is 25 %, the total compressibility 5.10- bar 1 and the fluid 
viscosity 1 cP. A production well BHP-PROl, with a radius of 7.85 cm, free of any skin 
effect, is at the centre of the reservoir : it is surrounded by four observation wells (BHP- 
OBS1, BHP-OBS2, BHP-OBS3, BHP-OBS4) (Figures 4). A numerical flow simulation 
allows to obtain, for this reservoir, a set of reference data comprising the bottomhole 
pressures of the five wells (Figure 5). 

The object of the inverse problem is to determine a reservoir model coherent with 
the pressure data, the permeabilities distribution being assumed to be unknown. Four 
optimization processes are therefore launched, starting from the same initial realization 
(Figure 4B). For each process, we consider a single optimization parameter, i.e. the 
deformation parameter. The first process (GDM1) takes up the conventional gradual 
deformation scheme with construction of a realization chain using the initial realization 
and a single complementary realization. The other three processes (GDMGBC3, 
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GDMGBC10 and GDMGBC30) illustrate the application of the method according to 
the invention : the chains are in this case elaborated with the initial realization and a 
composite realization obtained from the combination of 3, 10 and 30 complementary 
realizations. The composite realization is constructed as explained above (Eq.2) by 
integrating the information provided by the gradients of the objective function with 
respect to the Gaussian white noise. For each process, the evolution of the objective 
function (OF) as a function of the number k of flow simulations performed is shown in 
Figure 6. 

It can be observed that using the gradients and increasing the number of 
complementary realizations allows the objective function to decrease more rapidly for 
the same number of flow simulations. 



